#include <cmath>
#include <cstdlib>
#include <cstring>
#include <fstream>
#include <iostream>
#include <list>
#include <map>
#include <queue>
#include <string>
#include <vector>


#include <tbb/atomic.h>
#include <tbb/mutex.h>
#include <tbb/parallel_for.h>
#include <tbb/parallel_reduce.h>
#include <tbb/parallel_sort.h>
#include <tbb/task_scheduler_init.h>


typedef unsigned char Letter;

typedef std::string Sequence_name;
typedef std::string Sequence;
typedef std::vector<Letter> DNA_sequence;
typedef std::pair<Sequence_name, Sequence> Named_sequence;
typedef std::pair<Sequence_name, DNA_sequence> Named_DNA_sequence;

typedef std::list<Named_DNA_sequence> Input_sequences;

typedef std::vector<int> blocked_seq_t;
typedef std::pair<std::pair<size_t, size_t>, std::pair<size_t, size_t> > result_t;
typedef std::vector<result_t> results_t;
typedef std::pair<std::string, results_t> named_results_t;


struct DNA
{
	static Letter const A;
	static Letter const T;
	static Letter const G;
	static Letter const C;
	static Letter const NO_SYMBOL;
};

Letter const DNA::A = 0;
Letter const DNA::T = 1;
Letter const DNA::G = 2;
Letter const DNA::C = 3;
Letter const DNA::NO_SYMBOL = 4;

Letter charToDNA(char c)
{
	if (c == 'A')
	{
		return DNA::A;
	}
	if (c == 'T')
	{
		return DNA::T;
	}
	if (c == 'G')
	{
		return DNA::G;
	}
	if (c == 'C')
	{
		return DNA::C;
	}
	return DNA::NO_SYMBOL;
	return 0;
}


/* Represents block of _size characters. */
class block_t
{
public:
	static size_t const DEFAULT_SIZE;

	block_t(); // Default size
	block_t(size_t size);
	block_t(size_t size, size_t value);

	size_t value() const;
	void nextSymbol(Letter symbol);

	size_t maxSize() const;
	size_t size() const;

private:
	void setMasks();

private:
	size_t _size;
	size_t _current_size;
	size_t _value;

	size_t _extract_value_mask;
	size_t _change_value_mask;

	static size_t const LETTER_MASK;
	static size_t const SIZE_OF_LETTER;
	static size_t const MAX_SYMBOL_COUNT;
};

size_t const block_t::DEFAULT_SIZE = 10;
size_t const block_t::LETTER_MASK = 3;
size_t const block_t::SIZE_OF_LETTER = 2;
size_t const block_t::MAX_SYMBOL_COUNT = 4;

block_t::block_t()
	: _size(10)
	, _current_size(0)
	, _value(0)
{
	setMasks();
}

block_t::block_t(size_t size)
	: _size(size)
	, _current_size(0)
	, _value(0)
{
	setMasks();
}

block_t::block_t(size_t size, size_t value)
	: _size(size)
	, _current_size(size)
	, _value(value)
{
	setMasks();
}

size_t block_t::value() const
{
	return _value;
}

void block_t::nextSymbol(Letter symbol)
{	
	if (symbol < DNA::NO_SYMBOL)
	{
		if (_current_size < _size) 
		{
			++_current_size;
		} 
		_value <<= SIZE_OF_LETTER;
		_value += symbol;
		_value &= _change_value_mask;
	}
}

size_t block_t::size() const
{
	return _current_size;
}

size_t block_t::maxSize() const
{
	return _size;
}

void block_t::setMasks()
{
	_extract_value_mask = 1;
	for (size_t i = 0; i != SIZE_OF_LETTER * _size; ++i)
	{
		_extract_value_mask <<= 1;
		++_extract_value_mask;
	}
	_change_value_mask = 1;
	for (size_t i = 0; i != SIZE_OF_LETTER * (_size - 1); ++i)
	{
		_change_value_mask <<= 1;
		++_change_value_mask;
	}
}


Named_DNA_sequence & get_dna_sequence(std::ifstream & stream, Named_DNA_sequence & result_buff)
{
	result_buff.first.clear();
	result_buff.second.clear();
	Sequence_name & name = result_buff.first;
	DNA_sequence & sequence = result_buff.second;

	if (stream.get() == '>')
	{
		stream >> name;
	}

	while (!stream.eof())
	{
		int c = stream.get();
		if (c == '>')
		{
			stream.unget();
			break;
		}
		else
		{
			if (c == 'A')
				sequence.push_back(DNA::A);
			else if (c == 'T')
				sequence.push_back(DNA::T);
			else if (c == 'G')
				sequence.push_back(DNA::G);
			else if (c == 'C')
				sequence.push_back(DNA::C);
		}
	}

	return result_buff;
}





// All program logic is in processor_t class
class processor_t
{
public:
	processor_t(size_t threads_num, size_t min_match_len);

	void sortResults();
	std::map<std::string, results_t> const & results() const;

	void setReferenceString(DNA_sequence * seq);
	void setNewInputString(Named_DNA_sequence const & seq);

private:

	size_t _threads_num;
	size_t _min_match_len;
	size_t _block_size;

	std::map<std::string, results_t> _results;
	std::string _last_seq_name;


	DNA_sequence const * _reference_str;
	std::vector<size_t> _reference_blocked_str;
	std::vector<std::vector<size_t> > _reference_hash;



	results_t & processBlockedSequence(blocked_seq_t const & sequence, size_t seq_id
		, DNA_sequence const & seq, results_t & results_buff, tbb::mutex & m) const;

	void makeBlocks(std::vector<Letter> const & input, std::vector<size_t> & result);

	class reference_processing_func;
	class input_processing_func;
};

// Show results
std::ostream & operator<< (std::ostream & out, processor_t const & processor);


processor_t::processor_t(size_t threads_num, size_t min_match_len)
	: _threads_num(threads_num)
	,_min_match_len(min_match_len)
{
	size_t max_block = (_threads_num > 12) ? 8 : 6;
	_block_size = (min_match_len < max_block) ? min_match_len : max_block;
}

std::map<std::string, results_t> const & processor_t::results() const
{
	return _results;
}


class processor_t::reference_processing_func 
{
public:
	reference_processing_func(processor_t * processor, std::vector<tbb::mutex> * locks) 
		: _processor(processor)
		, _locks(locks)
	{}

	void operator() (tbb::blocked_range<size_t> const & r) const
	{
		block_t block(_processor->_block_size);
		size_t start = _processor->_block_size - 1 + r.begin();
		for (size_t i = r.begin(); i < start; ++i)
		{
			block.nextSymbol(_processor->_reference_str->at(i));
		}

		size_t end = r.end() + _processor->_block_size - 1;
		for (size_t i = start; i < end; ++i)
		{
			block.nextSymbol(_processor->_reference_str->at(i));

			_locks->at(block.value() % _locks->size()).lock();
			_processor->_reference_hash[block.value()].push_back(i - _processor->_block_size + 1);
			_locks->at(block.value() % _locks->size()).unlock();
		}
	}

private:
	processor_t * _processor;
	std::vector<tbb::mutex> * _locks;
};


void processor_t::setReferenceString(DNA_sequence * seq)
{
	_reference_str = seq;

	static size_t locks_count = _threads_num * 1024;
	size_t blocks_count = seq->size() - _block_size + 1;
	size_t size = 1 << (_block_size * 2);
	_reference_hash.resize(size);
	std::vector<tbb::mutex> locks(locks_count, tbb::mutex());

	size_t interval = std::ceil((double)blocks_count / (double)_threads_num);

	tbb::parallel_for(tbb::blocked_range<size_t>(0, blocks_count, interval), reference_processing_func(this, &locks));
}

class processor_t::input_processing_func
{

public:

	input_processing_func(processor_t const * processor, tbb::mutex * mutex, DNA_sequence const & seq, results_t * res)
		: _processor(processor)
		, _mutex(mutex)
		, _seq(seq)
		, _res(res) 
	{}

	void operator() (tbb::blocked_range<size_t> const & r) const
	{
		block_t block(_processor->_block_size);
		size_t start = _processor->_block_size - 1 + r.begin();

		for (size_t i = r.begin(); i < start; ++i)
		{
			block.nextSymbol(_seq[i]);
		}

		size_t end = r.end() + _processor->_block_size - 1;
		for (size_t i = start; i < end; ++i)
		{
			block.nextSymbol(_seq[i]);
			size_t idx_seq_start = i - _processor->_block_size + 1;

			if (!_processor->_reference_hash[block.value()].empty())
			{
				std::vector<size_t> const & lst = _processor->_reference_hash[block.value()];
				for (std::vector<size_t>::const_iterator it = lst.begin(); it != lst.end(); ++it)
				{
					if (idx_seq_start != 0 && (*it) != 0 && _seq[idx_seq_start - 1] == _processor->_reference_str->at((*it) - 1))
					{
						continue;
					}

					size_t idx_ref = (*it);
					size_t idx_seq = idx_seq_start;

					while (idx_ref < _processor->_reference_str->size()
						&& idx_seq < _seq.size()
						&& _processor->_reference_str->at(idx_ref) == _seq[idx_seq])
					{
						++idx_seq;
						++idx_ref;
					}

					--idx_seq;
					--idx_ref;

					if ((idx_seq - idx_seq_start + 1) >= _processor->_min_match_len && (idx_ref - (*it) + 1) >= _processor->_min_match_len)
					{
						_mutex->lock();
						_res->push_back(std::make_pair(std::make_pair(*it + 1, idx_ref + 1), std::make_pair(idx_seq_start + 1, idx_seq + 1)));
						_mutex->unlock();
					}

				}
			}
		}
	}

private:

	processor_t const * _processor;
	tbb::mutex * _mutex;
	DNA_sequence const & _seq;
	results_t * _res;
};


void processor_t::setNewInputString(Named_DNA_sequence const & seq) 
{
	_last_seq_name = seq.first;
	results_t & results = _results[_last_seq_name];

	size_t blocks_count = seq.second.size() - _block_size + 1;
	tbb::mutex m;

	tbb::parallel_for(
		tbb::blocked_range<size_t>(0, blocks_count),
		input_processing_func(this, &m, seq.second, &results)
	);
}

std::ostream & operator<<(std::ostream & out, processor_t const & processor)
{
	///+ show processor.result()
	for (std::map<std::string, results_t>::const_iterator it = processor.results().begin(); it != processor.results().end(); ++it)
	{
		out << it->first << "\n";
		results_t const & res = it->second;
		for (size_t i = 0; i < res.size(); ++i)
		{
			bool flag = true;
			size_t j = i + 1;
			while (flag && j < res.size()) 
			{
				if ((res[j].first.second - 1 == res[i].first.second 
					&& res[j].first.first - 1 == res[i].first.first
					&& res[j].second == res[i].second) 
					||
					(res[j].second.second - 1 == res[i].second.second 
					&& res[j].second.first - 1 == res[i].second.first 
					&& res[j].first == res[i].first))
				{
					flag = false;
				} 
				else if (res[j].first.second - 1 > res[i].first.second) 
				{
					break;
				}
				++j;
			} 

			if (flag) 
			{
				out << res[i].first.first << " " << res[i].first.second << " " 
					<< res[i].second.first << " " << res[i].second.second << "\n";
			}
		}
	}
	return out;
}

class results_order_t
{
public:
	bool operator() (result_t const & left, result_t const & right) const 
	{
		if (left.first.second != right.first.second)
		{
			return left.first.second < right.first.second;
		}
		return left.second.second < right.second.second;
	}
};

void processor_t::sortResults()
{
	for (std::map<std::string, results_t>::iterator it = _results.begin(); it != _results.end(); ++it)
	{
		tbb::parallel_sort(it->second.begin(), it->second.end(), results_order_t());
	}
}



void defaultArgs(int & argc, char ** & argv)
{
	argc = 5;
	argv = new char* [5];

	argv[1] = new char[3];
	strcpy(argv[1], "1");

	argv[2] = new char[3];
	strcpy(argv[2], "20");

	argv[3] = new char[10];
	strcpy(argv[3], "refseq.txt");

	argv[4] = new char[9];
	strcpy(argv[4], "input.txt");
}


class reader_t
{
	Named_DNA_sequence * _reference;
	std::vector<Input_sequences> * _sequences;
	char ** _files;
public:
	reader_t(Named_DNA_sequence * reference, std::vector<Input_sequences> * sequences, char ** files)
		: _reference(reference)
		, _sequences(sequences)
		, _files(files)
	{}

	void operator ()(tbb::blocked_range<size_t> const & r) const
	{
		for (size_t i = r.begin(); i < r.end(); ++i)
		{
			std::ifstream seq_file(*(_files + i), std::ios::in);
			if (i == 0)
			{
				get_dna_sequence(seq_file, *_reference);
			}
			else
			{
				Input_sequences & result = _sequences->at(i - 1);
				while (!seq_file.eof())
				{
					result.push_back(Named_DNA_sequence());
					get_dna_sequence(seq_file, result.back());
				}
			}
			seq_file.close();
		}
	}
};


int main(int argc, char ** argv)
{
	if (argc < 5) 
	{
		std::cout << "Incorrect arguments number" << std::endl;
		defaultArgs(argc, argv);
	}

	size_t worker_threads_num = atol(argv[1]);
	processor_t processor(worker_threads_num, atol(argv[2]));

	Named_DNA_sequence reference_DNA;
	std::vector<Input_sequences> input_DNA_sequences(argc - 4);

	tbb::task_scheduler_init init(std::min(worker_threads_num, (size_t)argc - 3));
	tbb::parallel_for(tbb::blocked_range<size_t>(0, argc - 3, 1), reader_t(&reference_DNA, &input_DNA_sequences, argv + 3));
	init.terminate();

	init.initialize(worker_threads_num);
	processor.setReferenceString(&reference_DNA.second);
	for (std::vector<Input_sequences>::const_iterator it = input_DNA_sequences.begin(); it != input_DNA_sequences.end(); ++it)
	{
		for (Input_sequences::const_iterator seq_it = it->begin(); seq_it != it->end(); ++seq_it)
		{
			processor.setNewInputString(*seq_it);
		}
	}

	processor.sortResults();
	std::cout << processor;

	return 0;
}
